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Abstract 

Several aspects of the theory of epitaxial crystal growth from atomic or molecular 
beams are developed from the perspective of statistical physics. Lectures are devoted 
to the rate equation theory of two-dimensional nucleation and its limitations; the 
growth of multilayer wedding cakes in the presence of strong step edge barriers; the 
continuum theory of mound coarsening; and growth-induced step meandering on 
vicinal surfaces. 
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1 Introduction 

The growth of a crystalline film from a molecular or atomic beam, commonly 
referred to as Molecular Beam Epitaxy (MBE), is a simple example of a self- 
assembly process. In contrast to crystallization from the melt, which often 
leads to dendrites and other ramified patterns [1], MBE growth can be de- 
scribed without reference to the transport of matter, latent heat or impurities 
in the fluid phase. The remarkable richness of patterns forming during MBE 
is determined solely by processes which occur locally at the surface. Moreover, 
in the case of homoepitaxial growth, in which a film is grown on a substrate 
of the same material, energetic determinants such as interfacial free energies 
and misfit strain are absent, so that the film morphology is governed primarily 
by growth kinetics. This makes homoepitaxy an ideal laboratory in which to 
study the emergence of mesoscopic patterns and their associated length scales 
through self-organization processes far from equilibrium. An additional benefit 
is that, since the invention of scanning tunneling microscopy, this laboratory 
is open to direct visual inspection. 
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The main steps in the growth process can be summarized as follows^. Atoms 
are deposited at rate F. They migrate along the surface with a two-dimensional 
diffusion coefficient D. When two atoms meet they form a dimer. Dimers may 
subsequently disintegrate, or they may grow by aggregation of further atoms 
into trimers and larger clusters. Once a substantial fraction of the surface is 
covered by two-dimensional island clusters, these begin to coalesce and a full 
atomic layer forms, on which the processes involved in producing the first layer 
repeat themselves. 
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Fig. 1. An atom descending from a step edge experiences an additional energy 
barrier A£g. 

At this stage of the growth process, the fate of atoms deposited on top of first 
layer islands becomes important. Such atoms may either descend from the 
island, thus contributing to the growth of the island edge, or they may remain 
on the island, promoting the nucleation of the next layer. On many crystal sur- 
faces, the diffusion of atoms between different atomic layers is suppressed due 
to additional energy barriers, which an atom crossing a step has to overcome. 
This phenomenon was first observed experimentally by Ehrlich and Hudda 
[2], and some of its consequences for the growth of stepped surfaces were an- 
alyzed by Schwoebel and Shipsey [3]. The atomistic origin of the additional 
step edge barrier is illustrated in Figure 1: An atom descending from a step 
edge passes through a transition state of very low coordination, which implies 
poor binding and thus a higher energy. This picture is oversimplified, because 
in many cases descent by concerted exchange is more facile than hopping [4]. 
Nevertheless, it is generally true that the rate for interlayer diffusion, in the 
following denoted by D', is smaller than the in-layer diffusion constant D. 

In-layer and interlayer diffusion, as well as all other atomic processes involved 

1 Throughout these notes, the unit of length will be the substrate lattice spacing. 
Thus the deposition rate F denotes the number of atoms deposited per unit time 
and adsorption site, and the diffusion coefficients D and D' are actually hopping 
rates. All three quantities have units of inverse time. 
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Fig. 2. Growth of Pt on Pt(lll) at T = 440K [5]. The total coverage is (a) 0.3 mono- 
layers (ML), (b) 3 ML, (c) 12 ML and (d) 90 ML. The image size is 2600A x 3450A. 
Courtesy of T. Michely. 

in the growth, decay and shape changes of two-dimensional islands, are ther- 
mally activated. Let us recall what this means, using the in-layer migration 
process as illustration. To a good approximation, the motion of an adsorbed 
atom (an adatom) on a crystal surface can be viewed as a two-dimensional 
random walk between adsorption sites. In hopping from one adsorption site 
to another, the adatom has to overcome an energy barrier Ed- The energy is 
provided by thermal substrate vibrations. This implies the familiar Arrhenius 
form 



for the diffusion coefficient. In (1), T is the substrate temperature, k B the 
Boltzmann constant, and D Q is an attempt frequency with a typical magni- 
tude around 10 13 s _1 . The main role of temperature in MBE growth is that 
it regulates, through expressions like (1), the relative rates of the different 
activated processes on the surface. 

In these lectures we will explore some of the pattern forming phenomena that 
arise through the interplay of the three kinetic rates D, D' and F . In Section 2, 
the classical atomistic rate equation theory of two-dimensional nucleation will 
be briefly reviewed. Its central result is a scaling law, Eq.(12), which relates 
the spacing between first layer island clusters to the ratio D/F. The limita- 
tions of the classical theory become evident in the treatment of nucleation 
on top of islands, which requires statistical arguments beyond rate equations. 
Sections 3 and 4 are devoted to the mound patterns which appear generi- 
cally in multilayer growth (see Figure 2 for an example). The approach taken 
in Section 3 is quite atomistic, focusing on the distribution of the deposited 
mass among the different layers and the corresponding mound shapes. Section 
4 provides a more macroscopic perspective, which allows one to address also 
the mass transport between mounds and the coarsening of the pattern. The 
final Section 5 is devoted to growth on vicinal surfaces, which are intrinsically 
anisotropic due to the presence of preexisting steps. 



D = D e~ ED/kBT 
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The treatment of these topics is far from exhaustive. More details can be 
found in several recent books and review articles [6-12]. The lectures do not 
cover the theory of kinetic roughening by stochastic fluctuations, which in 
fact initiated, through the seminal work of Kardar, Parisi and Zhang [13], 
the interest of the statistical physics community in film growth. The reasons 
for this omission are twofold. First, the subject has been extensively reviewed 
elsewhere [11,12,14,15]. Second, despite its tremendous conceptual impact, the 
relevance of kinetic roughening theory to real film growth has still not been 
quantitatively demonstrated^]. This is in contrast to the growth instabilities 
discussed in the present lectures, where a quantitative linking of mesoscopic 
patterns to specific atomistic processes appears to be well within reach. A 
possible manifestation of kinetic roughening in a setting which is relevant to 
MBE growth is the noise-induced damping of growth oscillations in layer-by- 
layer growth. Here kinetic roughening theory has been used to establish scaling 
relationships between the damping time and the ratio D/F [11,17-19]. 



2 Two-dimensional nucleation 

The formation of an atomic layer on a high symmetry substrate without steps 
or defects has to proceed through the congregation of mobile adatoms into 
stable clusters, which subsequently grow by accretion. The process is analo- 
gous to the nucleation of a condensed phase out of a supersaturated gas, as 
described by thermodynamic nucleation theory [20]. The central object of this 
theory is the critical nucleus, which defines the free energy barrier that has 
to be surmounted to reach the stable phase. The size of the critical nucleus is 
inversely proportional to the supersaturation. 

In far from equilibrium growth the critical nucleus size may reach atomic 
dimensions, thus precluding a straightforward application of thermodynamic 
concepts and necessitating the development of an atomistic theory of nucle- 
ation kinetics [21]. This was achieved in the 1960's and 70's by Zinsmeister, 
Stowell, Venables and others; extensive reviews are available [22-24]. Precise 
experimental tests of the theory have become possible only recently [25,26]. 
In addition, large scale computer simulations play an increasingly important 
role in establishing the validity and limitations of nucleation theory, as they 
allow for separate scrutiny of the various assumptions going into the theory. 

The atomistic theory of two-dimensional nucleation is summarized in the next 
section. We then turn to the problem of second layer nucleation on top of 
islands, which differs from the nucleation of the first layer because of the 
confinement of the atoms by step edge barriers. This problem is of conceptual 

2 Experimental work until 1995 has been reviewed in [16]. 



4 



interest because it illustrates the limitations of mean field rate equation theory 
[27-29]. At the same time a quantitative theory of second layer nucleation is a 
necessary prerequisite for the discussion of multilayer growth in the subsequent 
lectures. 



2.1 Rate equation theory 



The classical approach to nucleation kinetics starts from balance or rate equa- 
tions for the areal concentrations n s of clusters consisting of s atoms; ri\ is the 
adatom density, 112 the density of dimers, and so on. To be precise, we define 
n s as the number of clusters per surface area, averaged over a region con- 
taining a large number of clusters. If the adatoms are the only mobile species 
(the mobility of larger clusters is negligibly small), then clusters grow solely 
by aggregation of single adatoms. Defining T s to be the net rate at which 
s + 1-clusters form from s-clusters, we have for s > 2 

^ = r s _!-r s ( a >2). (2) 



The net formation rates T s can be written as 

T s = cTsDrixris - 7 s+ in s+ i, (3) 



where 7 S is the rate at which adatoms detach from a s-cluster and the di- 
mensionless capture number a s accounts for the propensity of a s-cluster to 
absorb an adatom (see below). The chain (2) of aggregation equations is fed 
by the adatom density ri\. If desorption from the surface can be neglected (the 
complete condensation limit [22,24]), adatoms are lost only through dimer for- 
mation and capture at larger clusters, and the adatom rate equation reads 

^ = F- 2ri -£r, (4) 

0,1 s>2 



The deposition rate F is defined as the number of atoms arriving per unit 
time and surface area. 

In principle, Eqs. (2-4) provide a complete description of the nucleation process; 
in practice, they contain far too many (generally unknown) kinetic parame- 
ters to be useful. This difficulty is commonly circumvented by introducing a 
distinction between stable and unstable clusters, and postulating a separation 
of time scales between the kinetics of the two kinds. Stable clusters of sizes 
s > + l are assumed not to decay, i.e. 7 S = for s > i*, while the detachment 
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of adatoms from unstable clusters with s < i* occurs sufficiently rapidly to 
establish thermodynamic equilibrium between the different size populations. 
It is important to note that, in contrast to thermodynamic nucleation the- 
ory the critical cluster size i* introduced here contains a kinetic element, as 
it refers to stability and equilibration only on the time scale relevant to the 
deposition experiment. 

Next the total density N of stable clusters, also referred to as islands in the 
following, is introduced through 

oo 

N= £ n s . (5) 

s=i* + l 



Summing (2) from s = i* + I this is seen to evolve according to 

— = ai*Dnirii*. (6) 



The assumption of thermal equilibrium among unstable clusters implies that 
the net formation rates Y s vanish for 1 < s < i* — 1, and induces the detailed 
balance relations [21] 

n s ss n{e Ea/kBT (2 < s < i*) (7) 



between the unstable cluster concentrations and the adatom concentration. 
Here E s is the total (positive) binding energy of an s-cluster, i.e. the energy 
needed to disperse the cluster into single adatoms; note that E\ = 0. 

Using (7) the nucleation rate on the right hand side of (6) can be expressed in 
terms of the adatom density. To complete the description, the rate equation for 
rii is simplified by introducing the average capture number for stable islands 

oo 

a = N- 1 ]T n s a s . (8) 

s=i*+l 



Equation (4) can then be written in the form 

dni 



F — (j^Dnirii* — aDriiN. (9) 
dt 



Together (6), (7) and (9) form a closed set of equations from which the island 
and adatom densities can be computed. 
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Before turning to the solution of these equations, some remarks concerning the 
capture numbers o~ s are in order. In the early literature on atomistic nucleation 
theory, a geometric view of the capture numbers as effective cross sections for 
atoms colliding with clusters prevailed. More appropriately, they are defined 
through the requirement that the flux of adatoms to a cluster of size s should 
be, on average, equal to o s Dn\. The capture number of a cluster evidently 
depends on its size, but in general it depends on the sizes and locations of 
surrounding clusters as well, because these affect the adatom concentration 
field. The calculation of a s involves the solution of a diffusion equation for 
the adatom concentration, with appropriate boundary conditions represent- 
ing both the capture of adatoms at the cluster of interest, and the presence 
of other clusters far away [30] . In this sense the capture numbers retain some 
information about the spatial arrangement of the clusters, which is otherwise 
not accounted for in the rate equations; in the jargon of statistical mechan- 
ics, the rate equation description constitutes a mean field approximation. In 
practice, both o> and a turn out to be slowly varying functions which can 
be replaced by constants for many (though not all) purposes. We adopt this 
simplification for the present discussion. 

The solutions of the coupled equations (6) and (9) display two temporal 
regimes. In the early time, transient nucleation regime the loss terms on the 
right hand side of (9) are negligible. The adatom concentration increases pro- 
portional to the total coverage 6 = Ft, and the island density grows rapidly 
as N ~ O l * +2 . This regime ends when capture of adatoms at stable islands 
becomes appreciable, at a coverage which can be estimated by comparing the 
first and last terms on the right hand side of (9). In the subsequent steady 
state regime these two terms balance completely, and the adatom density is 
determined by the island density through 



Inserting this into (6,7) and integrating in time yields the central result of 
nucleation theory, 



The most important feature of this expression is that it takes the form of a 
scaling relation 






(12) 
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Fig. 3. Kinetic processes involved in second layer nucleation. The figure shows a 
circular island viewed from the side. 

between the island number density and the ratio of the two basic kinetic rates 
D and F of the deposition process, with the scaling exponent \ taking the 
value x — i* /{i* + 2) (for a generalized expression see (24)). 

The rate equations formulated above are restricted to the low coverage regime, 
since finite coverage effects such as direct impingement and cluster-cluster 
coalescence have been neglected. In the absence of cluster mobility, coalescence 
provides the only mechanism by which the island density can decrease. This 
produces a maximum in N(Q) which is often used as a convenient reference 
point for the experimental determination of cluster densities. 



2.2 The rate of second layer nucleation 



We now ask how the nucleation process is modified when it occurs on top of an 
island, rather than on the (unbounded) substrate. The geometry is illustrated 
in Figure 3. Atoms are being deposited at rate F onto a circular island of 
radius R. They diffuse on the island with an in-layer diffusion constant D, 
and descend from it with an (average) interlayer diffusion rate D' < D. We 
are looking for the probability per unit time, u, for a nucleation event to 
occur on the island. Assuming that dimers of adsorbed atoms are stable (as 
they will be at sufficiently low temperatures), nucleation takes place as soon 
as two adatoms meet. In the terminology of the preceding subsection, this 
corresponds to i* = 1; for a generalization of the following considerations to 
i* > 1 see [29]. 

The kinetic rates F, D and D' combine to form three relevant time scales: 
The interarrival time 
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between subsequent depositions onto the island; the diffusion time 

r D ~ R 2 /D (14) 



required for an atom to diffuse once across the island, or for two atoms to meet; 
and the residence time r which an atom spends on the island before descend- 
ing, when no nucleation event takes place. The calculation of the residence 
time requires the solution of a stationary diffusion problem with appropriate 
boundary conditions at the step edge, which yields [27,31] 

t= 1 . (15) 

8D 2D' K J 

The residence time is comparable to the diffusion time td if the suppression 
of interlayer transport is weak, in the sense that D / D' <C R, while in the 
opposite regime of strong step edge barriers, D/D' ^> R, the residence time 
becomes r = R/2D' independent of v. 

We will focus on the latter regime in the following. Then r ^> td, which implies 
that two atoms are certain to meet once they are present simultaneously on the 
island. The probability for this to occur is r/(r + At) t f At if At 3> r, which 
is true for reasonable deposition fluxes. Multiplying this by the total number 
of atoms deposited onto the island per unit time, we obtain the expression 
[27] 

uj= wr = ^- (16) 



for the nucleation rate, which is exact under the stated conditions. 



2. 3 The limitations of mean field nucleation theory 



The first calculation of the rate of second layer nucleation [31] was based on 
the rate equation (6), in which the average adatom density rt\ is replaced by 
its local value n(r) at a point f on the island. The local nucleation rate I(r), 
which counts the number of nucleation events per time and lattice site, is then 
given by 

/(f) = a i *De E >* /kBT n(rf +1 . (17) 



Computing the adatom density profile from the stationary diffusion equation 
and integrating over the island area, the total nucleation rate on a circular 
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island is obtained. Comparison with the exact expression (16) for %* — 1 shows 
that the rate equation approach overestimates (16) by a factor of the order of 
D/RD' > 1 [27]. 

To elucidate the origin of this discrepancy, it is useful to view the encounter 
of the two atoms as a first passage problem in the four-dimensional space 
spanned by their coordinates. The joint diffusion of the atoms stops either 
when the two reach adjacent lattice sites (nucleation) or when one of them 
escapes, whatever happens first. The detailed analysis of this process [32] 
reveals the nature of the approximation based on (17): Setting the nucleation 
rate proportional to the square of the adatom density amounts to treating the 
atoms as noninteracting, in the sense that they are allowed to continue their 
diffusion even after they have met; in this way a single pair can accumulate 
several (spurious) nucleation events, and the nucleation rate is overestimated. 

For a quantitative comparison, we note that (for %* = 1) the nucleation proba- 
bility p nuc per atom can be written quite generally as the product of the (mean) 
adatom density and the number N dis (r) of distinct sites the atom encounters 
during its lifetime [33]. The validity of this statement becomes evident by as- 
suming that all other adatoms are immobile: Then it is clear that repeated 
visits to the same unoccupied site do not increase the chance for nucleation. 
The lifetime of the atom equals its residence time in the case of second layer 
nucleation, but the argument applies equally well to the nucleation of the first 
layer, where the lifetime is determined by capture at stable islands, see below. 

Using the general relation [27] 

t = n/F (18) 

between the residence time and the mean adatom density n, which follows from 
simple mass balance considerations, the total nucleation rate on the island can 
then be written as 

U = (AtrVnuc = FN ^t 

In the strong barrier limit, iV dis equals the total number of sites on the island, 
and (19) reproduces (16). To compare (19) to the rate equation approximation, 
we multiply the local nucleation rate (17) by the island area and obtain, in 
order of magnitude, 

uj mi ~ R 2 D(n) 2 ~ FN^, (20) 

where (18) has been used, and the number 7V a n = Dr of all sites visited by an 
adatom during its residence has been introduced. Thus the expressions (19) 
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and (20) differ, in general, by a factor N a \\/N dis > 1. 

The distinction between AT all and affects also the mean field calculation 
of the density of first layer islands in Section 2.1 [33]. In this case the life time 
of a freshly deposited adatom is of the order r ~ 1/(ND), and the average 
adatom density is ri\ = Ft ~ F/ND, compare to (10). The theory of two- 
dimensional random walks provides the expression iVdis ~ 7iDr/ln(Dr) for 
Dt ^> 1 [34]. The actual value of N is fixed by the requirement that, out of 
the approximately 1/N adatoms deposited in the area occupied by one island, 
only one (or two!) forms a nucleus, i.e. that p nuc ~ niN dis (r) ~ N. This yields 
finally 

AT 3 ln(l/A0^, (21) 



which coincides with the scaling law (12) for %* — 1 only up to a logarithmic 
correction; the leading behavior of (21) for D/F ^> 1 is 

N ~ {F/D) 1/3 [\og{D/F)]- 1/3 . (22) 



In a power law fit of N versus Fj D this will tend to produce scaling exponents 
which are smaller than 1/3. 

The logarithmic factor can be reproduced within rate equation theory by using 
the expression 

4tt 1 

en ~ ~ (23) 

1 \n[(D/F)m] ln(iV) 1 ' 



for the capture number of adatoms [30]. Reduced rate equations of the form 
(6,9) which incorporate the logarithmic correction are able to quantitatively 
reproduce the results of kinetic Monte Carlo simulations [35]. 

The random walk picture has also been helpful in clarifying the case of one- 
dimensional diffusion. In one dimension the distance between islands is iV -1 , 
and correspondingly the adatom lifetime in the steady state regime is of 
the order of r ~ 1/(N 2 D) and the adatom density is m ~ F/(N 2 D). The 
number of distinct sites visited by a one-dimensional random walk grows as 
Ndi S (T) ~ V Dr. Inserting these expressions into the condition p nuc ~ N yields 
the estimate iV ~ {F/D) 1 ^ for i* — 1 [33]. Similar considerations can be 
employed to derive a useful general formula for the exponent x i n the scaling 
law (12) [29,36]. It reads 



d + 2i* + mm[di*,2]' 
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Here d = 1, 2 denotes the dimensionality of diffusion. Logarithmic corrections 
similar to (21) arise whenever di* = 2, because di* is the effective dimension- 
ality of the random walk in configuration space which describes the nucleation 
process, and two-dimensional random walks are well known to be marginally 
space filling [34]. 



3 Wedding cakes 

It was first predicted by Villain [37] that the existence of step edge barriers 
- the fact that D' < D - implies a growth instability, in which a mound 
morphology develops on a flat crystal surface. This phenomenon, which was 
subsequently observed on a variety of substrates [38-43], is the subject of the 
following two lectures. 

3.1 Poisson growth 

To appreciate the relevance of interlayer transport for multilayer growth, we 
first consider the case where it is completely absent, i.e. we set D' = 0. Then 
each adatom remains in the layer in which it was first deposited, and is in- 
corporated into that layer at an ascending step edge. This implies a simple 
evolution of the layer coverages 9 n , < 9 n < 1, where n = 1, 2, 3. ..counts the 
layers and 9 = 1 describes the substrate. The rate at which layer n grows is 
proportional to the exposed coverage 

<Pn-l = On-l - On (25) 

of the layer n — 1 below, and therefore the 9 n satisfy [44] 

^ = F(9 n _ 1 - 9 n ) (26) 

with the initial conditions 9 n >i(t = 0) = 0. It is straightforward to check that 
the solution reads 

n— 1 r\k 

9 n = l-e- e Yl ^, (27) 
k=o K - 

where 

oo 

= ]T 9 n = Ft (28) 

n=l 
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is the total deposited coverage. Correspondingly the exposed coverages follow 
a Poisson distribution with parameter 6, 

e -e@n 

y?n = — j— ■ (29) 



Since ip n is the probability that an arbitrary point on the surface resides on 
layer n, it can also be viewed as the probability distribution of the local film 
height, measured in units of the layer thickness. The mean height is 0, and the 
standard deviation of the <p n defines the surface width W, a common measure 
of film roughness, through 

oo 

w* = 5>-e)V (so) 

n=0 



For the Poisson distribution (29), the variance is equal to the mean, and 
therefore 

W = V& (31) 



for growth without interlayer transport. Remarkably, the expression (31) is 
independent of the in-layer diffusion rate. It represents the maximum surface 
roughness that can be generated by the randomness in the deposition flux, and 
is referred to in the literature as the statistical growth or random deposition 
limit [12]. Interlayer transport is solely responsible for reducing the roughness 
below this limit. 

While in-layer diffusion does not affect the vertical surface morphology, as en- 
coded in the layer coverages 9 n , it is certainly reflected in the lateral mass dis- 
tribution along the surface. This is illustrated in Figure 4 by a one-dimensional 
simulation. It shows the emergence of a fairly regular pattern of mound-like 
surface features with a characteristic pointed shape. Each mound consists of 
a tapering stack of islands upon islands, reminiscent of a wedding cake. In the 
following we argue that many properties of this pattern follow immediately 
from the expression (27) for the layer coverages [45] . 

Figure 4 clearly demonstrates that the lateral positions and sizes of the mounds 
originate in the growth of the first layer: The wedding cakes grow on the 
templates of the first layer islands, and their spacing is simply determined by 
the density N of first layer nuclei, which in one dimension is proportional to 
(F/D) 1 ^ 4 (see Sect. 2. 3). The persistence of the first layer pattern throughout 
the deposition of hundreds and thousands of layers requires, first, that no 
additional mounds nucleate during the later stages of growth, and, second, 
that neighboring mounds do not merge. The first requirement is met because 
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Fig. 4. Surface morphology in a one-dimensional growth model without interlayer 
transport, (a) Morphology after deposition of 1, 5.6, 16 and 32 monolayers, (b) 
Morphology after deposition of 45.25, 90.5, 181, 256 and 362 ML. The ratio of 
diffusion to deposition rate is D/F = 5 x 10 6 [45] 

the step spacing on the sides of the mounds decreases with increasing coverage 
(the mounds steepen), and therefore nucleation on the vicinal terraces making 
up the sides becomes highly unlikely. The merging of mounds is suppressed 
because the lateral positions of the maxima and minima of the surface profile 
are strongly correlated from one layer to the next. For the maxima this reflects 
the fact that nucleation occurs only in the top layer of each mound, while 
for the minima it can be attributed to the Zeno effect, the appearance of 
deep crevices between neighboring mounds [46]. Such a crevice closes very 
slowly, because fewer and fewer deposited atoms find their way to its bottom 
terrace. In particular, according to (27) the exposed fraction of the substrate 
ip — 1 — = e _e approaches zero only asymptotically, but does not vanish 
at any finite time. 

The resulting picture of the pattern forming process is readily generalized to 
growth on real, two-dimensional substrates. The nucleation of the first layer 
islands partitions the substrate into capture zones. Each zone supports a single 
mound, which is fed by the atoms deposited into the zone. This implies that 
the typical shape of the mounds can be read off from the layer distribution 
(27): The area A n of the n'th layer of a mound will be equal to 9 n A , where A 
is the area of the corresponding capture zone. A more transparent form of (27) 
is obtained in the limit of thick layers, ^> 1, when the Poisson distribution 
(29) can be replaced by a Gaussian of width v 7 ©. The layer distribution then 
follows by integration, 

e n = Q((n-G)/Ve) (32) 
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where 



*(*) = 1 - i / d V e ' y2/2 = \\ l ~ erf(*/v^)] (33) 

" - -oo 



and erf(s) denotes the error function. Equation (32) shows that the mounds 
attain a time-independent limiting shape when rescaled vertically by W = 

Ve. 



3.2 An improved model 



The main features of the Poisson growth model agree well with deposition 
experiments on Pt(lll) carried out for film thicknesses ranging from 0.3 to 
300 ML [5] (Figure 2). The mound shape at 130 ML, obtained by averaging 
the layer coverages over several mounds, matches the predicted shape func- 
tion (33) both in the valleys and on the slopes of the mounds, but differs 
significantly near the tops: Instead of the pointed peaks seen in Figure 4, the 
real mounds terminate in flat terraces of a characteristic lateral size (Figure 
5). This discrepancy should be no surprise, since the model assumption of 
zero interlayer transport becomes extremely unrealistic for an adatom that 
has been deposited onto the freshly nucleated, small top island of a mound. 
In the absence of ascending step edges, such an atom has no choice but to 
interrogate the island edge many times and eventually cross it, even if the 
crossing probability in each attempt is very small. As the island grows, the 
residence time of the adatom increases, until the probability for two atoms 
to be present simultaneously on the island becomes appreciable, and a new 
nucleus is formed on top of the island. 

We now describe a simple modification of the Poisson model which takes into 
account the delayed nucleation of the top layer [47] . A mound is approximated 
as a stack of concentric, circular islands. The base is an island of radius i?o 
which does not grow, and the radius of the n'th island is R n = Since 
there are no sinks for atoms on the top layer, all atoms deposited there have 
to attach to the descending step. The top terrace therefore absorbs all atoms 
landing on layers n top and n top — 1, and grows according to 



d0. 



it = F ^ top - 1 ' (34) 



while (26) still applies for n < n top . Equations (26,34) have to be supplemented 
by a rule for the nucleation of a new top terrace. A simple choice would 
be to posit that nucleation occurs whenever the current top layer reaches 
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Fig. 5. Experimentally determined mean mound shape (full line) compared to the 
predictions of the Poisson growth model (dashed line), and a fit to the shape func- 
tion in the improved model which includes nucleation on the top terrace (dotted 
line). The figure shows a mound composed of circular islands seen from the side. 
The experimental data are taken from [5], and are courtesy of T. Michely. The 
experimental conditions correspond to those of Figure 2 after deposition of 130 
monolayers. 

some critical coverage. More realistically, nucleation should be treated as a 
stochastic process governed by the nucleation rate oo computed in Section 
2.2. A deterministic rule which is close in spirit to stochastic nucleation can 
be obtained as follows. Suppose the current top terrace has nucleated at time 
t — 0, and denote its radius by i? top (t). Then the probability that no nucleation 
has occurred on the top terrace up to time t is given by [27] 



and the probability density of the time t of the next nucleation event is 
—dPo/dt. It follows that the mean value of Pq at the time of nucleation satisfies 



and thus Pq = 1/2. In the numerical implementation, we therefore monitor 
the increase of P during the growth of the top terrace and create a new top 
terrace when P = 1/2. For the nucleation rate u the expression (16) will be 
used. 

A full analytic solution of the model is difficult because the nucleation proba- 




(35) 




o 



(36) 
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Fig. 6. Convergence of numerically generated wedding cakes to the asymptotic 
shape. 



bility (35) depends on the size of the top terrace, which is determined by the 
entire history of the wedding cake through the coupled equations (26). The 
numerical solution of the model equations shows that, like in the case of Pois- 
son growth, the height profile of the mound converges to a time-independent 
asymptotic shape, when viewed relative to the mean film height Ft and 
rescaled horizontally by \J~F~i (Figure 6). For large coverages this implies a 
separation of time scales between the mound sides, which evolve slowly at 
typical step velocities of order l/y/Ft, and the top terrace, which reaches a 
coverage of order unity during the growth of one layer. It is therefore reason- 
able, as a first approximation, to assume that the former top terrace ceases 
entirely to grow once a new island has nucleated on top of it. This should 
produce a lower bound on the island radii reached at a given time. 

Denote by R^-i the radius of island n — 1 at nucleation of island n, and by t n 
the time of this nucleation event. Then during its tenure as the top terrace, 
the radius of island n grows, according to (34), as 

Rn(t) = V^(*-*n)i#fl. (37) 



The time t n +i of the next nucleation event is determined by evaluating (35) 
using (16) and (37), and setting the result equal to 1/2, 



Po(t n +i) = exp 



! (C- P i) f 



2D' 



J dt [F{t - t n )f 2 = 1/2. (38) 
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This implies a recursion relation 

< op = r 5 J 7 (rT,) 2/7 (39) 



for the size of the top terrace at nucleation of the next layer, where we have 
introduced the characteristic radius 



The time interval r n = t n+ i — t n during which n top = n is related to R l ° p 
through (37), and correspondingly satisfies 

r n+1 = F-\Fr n f/\ (41) 



It is easy to check that the recursion relations (39) and (41) approach fixed 
point values 

iC P - Re Fr n - 1 (42) 



exponentially fast in n. Thus asymptotically there is one nucleation event 
during the growth of one monolayer, and nucleation occurs when the radius of 
the top terrace has reached the value R c . The numerical solution shows that 
these statements remain valid for the full dynamics, although the approach of 
R l ° p and r n to their asymptotic values is slower than exponential due to the 
coupling to the lower layers (the deviations decay as 1/VFt). 

To derive the asymptotic mound shape analytically, we insert the ansatz (32) 
into (26) and expand for large 0. We find that the shape function $(s) has to 
satisfy the differential equation 

$"( s ) = - s $'( s ). (43) 



This shows that the inflection point of the profile, where $" = 0, is always 
located at s — 0, i.e. at n — Ft. The solution of (43) which satisfies the 
boundary condition lim s ^_ 00 $(s) = 1 reads 

$(s) = l-C[l + erf(s/>/2)], (44) 



where C is a constant of integration. The profile is cut off at the rescaled height 
Smax of the top terrace, where the coverage takes the value 9 C = (R c /R ) 2 , 

$(s max ) = Be- (45) 
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Accordingly, the height of the top terrace above the mean film thickness Ft is 
■SmaxV^f^- The Poisson shape (33) is recovered in the limit 9 C — > 0, s max — > oo, 
C -> 1/2. 

The two parameters C and s max of the shape function are related by (45), 
but to fix both of them, a further relation is required. This is provided by the 
normalization condition (28) which, using (32), translates into 

5 max 

j ds = J ds (1 - $(s)). (46) 

-oo 

Together equations (44,45,46) define a family of shape functions parametrized 
by 9 C . Other features of the mound shape can be easily extracted. For example, 
the surface width (30) turns out to be given by 

W 2 « (1 - 9 c )Ft (47) 

asymptotically. By fitting the experimental mound shapes to the predicted 
shape function, the top terrace size R c and hence, through (40), the interlayer 
diffusion rate D' can be determined [27]. Assuming equal preexponential fac- 
tors for in-layer and interlayer diffusion, the fit shown in Figure 5 yields an 
additional step edge barrier of AE S « 0.21 eV. 



3.3 The growth-induced current 

The approach of the preceding subsections has provided us with a fairly ac- 
curate description of the shape of individual mounds. However, it does not 
capture global features of the morphology, such as the spatial organization of 
the mounds and their size distribution. Most importantly, in many (though not 
all [5]) experiments the mounds are observed to coarsen, i.e. their typical lat- 
eral extent increases with film thickness [40-43]. This requires mass transport 
between mounds, and hence the basic assumption of the wedding cake model - 
that the entire mass deposited within the capture zone of a mound contributes 
to its growth - cannot be uphelcQ So far the only theoretical approach for 
treating coarsening is based on phenomenological continuum equations, which 
will be discussed in detail in the next section. Here our purpose is to lend 
some credibility to this approach by showing how the asymptotic shape of the 
wedding cakes can be derived within the phenomenological framework. 



See however Sect. 4. 6 for a mechanism of mound coarsening without lateral mass 
transport. 
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Fig. 7. Illustration of the uphill surface current generated by step edge barriers. 

The continuum theory of mound formation is based on the notion of a growth- 
induced, uphill surface current [6,11,37,48]. Figure 7 illustrates the basic idea: 
Because atoms deposited onto a vicinal terrace^] on the side of a mound attach 
preferentially to the ascending step, they travel on average in the direction of 
the uphill slope between the point of deposition and the point of incorporation. 
It is important to realize that this does not require atoms to actually climb 
uphill across steps [11,48]. 

The evaluation of the current is particularly simple when interlayer transport is 
completely suppressed, and nucleation on the vicinal terrace can be neglected. 
An atom deposited onto a vicinal terrace of width I then travels an average 
distance 1/2 to the ascending step, and hence the current is Fl/2. Describing 
the surface profile by a continuous height function h(f, t) which measures the 
film thickness (in units of the monolayer thickness) above a substrate point r, 
the local terrace width is I = |V/i| _1 , and hence the current is given by the 
expression 



The surface profile evolves according to the continuity equation 

dh 

-g + V • J = F. (49) 



Here we are specifically interested in radially symmetric mounds. Rewriting 
(49) in polar coordinates and using (48) yields the following evolution equation 
for the mound profile h(r,t): 



dh 
dt 




(50) 



4 A vicinal terrace is a terrace which is bounded by an ascending step on one side 
and a descending step on the other. For further discussion of vicinal surfaces see 
Section 5. 
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We want to show that (50) possesses solutions corresponding to the asymptotic 
mound shape derived in Section 3.2. To this end we make the ansatz 

h(r,t) = VFtip(r) + Ft. (51) 
Inserting this into (50) yields the differential equation 

i fe) = ^ (52) 

In order to establish the equivalence between the shapes described by ip(r) 
and the scaled coverage distribution $(s) of Section 3.2, we need to verify 
that the function ip(r) defined implicitly by 

(^) 2 = $W (53) 
solves (52). Indeed, taking the derivative with respect to tp on both sides yields 

r [d^y = _cRi e _, V2 (54) 

which reduces to (52) upon taking another derivative with respect to r. 

This calculation shows explicitly that the sides of the mounds evolve accord- 
ing to the continuum equation (49); in [45] the same result was obtained for 
a one-dimensional geometry. The continuum description does however not in- 
clude the nucleation on the top terraces, which has to be added to the profile 
determined by (52) as a boundary condition. The development of a continuum 
theory of epitaxial growth which explicitly incorporates nucleation remains a 
challenge for the future. 



4 Phenomenological continuum theory of mound formation 

4-1 Motivation of the evolution equation 

As it stands, the continuum equation (49) cannot be used globally, because 
the surface current (48) diverges for |V/i| — > 0. For small slopes the nucleation 
of islands on the vicinal terraces can no longer be neglected. In quantitative 
terms, island nucleation sets in when the width I of the vicinal terrace becomes 
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comparable to the distance Id between islands nucleated on a flat surface. 
Using (12), the island distance is estimated as 

l D « N- 1/2 ~ (D/F) x/2 , (55) 



with x — 1/3 for i* = 1. The nucleated islands capture part of the deposited 
atoms, thus reducing the flux to the ascending step, and hence the uphill 
current. It is easy to show that the resulting net current vanishes linearly for 
|V/i] — > [6,11]. Ignoring effects of crystal anisotropy (see Section 4.4), the 
current always points in the direction of V/i. It is then generally of the form 

j(V/i) = f(\Vh\ 2 )Vh, (56) 



where the function f(u) approaches a constant for u — > 0. To match the form 
(48) for large slope, we further require that / ~ u~ l for u — > oo. The transition 
between the two regimes should occur at a slope of the order of 1/Id- Under 
suitable rescaling [49] this slope, as well as all other dimensionful parameters 
can be set to unity, and the overall behavior of the current can be represented 
by the interpolation formula [39] 

/!(«) = -L_ (model I). (57) 
1 + u 



For large slopes the current functions (48) and (57) decrease indefinitely, but 
the uphill current remains nonzero for all slopes. This implies that mass is 
continuosly transported uphill, leading to unbounded steepening of the mor- 
phology. Experimentally, it is often observed that the mound slopes approach a 
constant, "magic" value after a transient phase of steepening, a phenomenon 
known as slope selection. This can be incorporated into the continuum de- 
scription by letting the current function vanish at some nonzero slope [50]. 
Scaling the selected slope to unity, a simple choice of the function / in (56) 
that incorporates slope selection is 

/ n (u) = l-« (model II). (58) 



The growth equations with current functions (57) and (58) will be referred to 
as model I and model II in the following. 

The model defined by Eqs. (49,56) still does not constitute a useful description 
of the growing surface. To understand what difficulty remains, let us linearize 
(49) around the flat solution, writing 

h(f,t) = Ft + e(r,t). (59) 
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We find a diffusion equation for e, but with a negative diffusion coefficient 
— /(0) [37]. This confirms that the flat surface is unstable, but the instabil- 
ity becomes arbitrarily violent on arbitrarily short length scales. To cure this 
unphysical feature, we have to introduce a smoothening term which counter- 
acts the uphill current on short length scales. Near thermal equilibrium, such 
smoothening is provided by the well known Gibbs-Thomson effect, which in- 
duces mass transport from positively curved regions of the surface (hilltops) to 
negatively curved regions (valleys). To leading order in a gradient expansion, 
this implies a mass currrent 

/smooth = KV{V 2 h), (60) 



where K is proportional to the product of the surface free energy and the 
adatom mobility [51,52]. The Gibbs-Thomson effect may not be directly rele- 
vant under the far-from-equilibrium conditions of MBE, but other mechanisms 
have been suggested that give rise to a smoothening current of the same gen- 
eral form [6,53,54]. 

We follow the common practice and add the divergence of Eq.(60) to the right 
hand side of (49). Repeating the linear stability analysis, we find exponentially 
growing or decaying perturbations e(r, t) ~ exp[i(q • r) + cr(q)t}, where the 
growth rate a of a perturbation of wavevector q is given by 

a(q) = f(0)\q\ 2 -m 4 - (61) 



The instability is now limited to an unstable band < \q\ < q c = y f(0)/K. 
Wavelengths below 27r/g c are stabilized by the smoothening current (60). The 
most rapidly growing mode, which dominates the initial pattern that emerges 
from a generic random perturbation, has wavelength^ 27r\/2/g c and a finite 
growth rate a(q c /\^2). The unphysical features of the instability have thus 
been successfully removed. 

As before, the coefficient K in (60) can be set to unity by suitable rescaling. 
In addition, the constant deposition rate F on the right hand side of (49) can 
be removed by letting h — > h — Ft. The final outcome of these considerations 
is the evolution equation 

|_ = -V-[/(|V/f)V/>]-(V 2 ) 2 /> (62) 
which will be analyzed in the remainder of this section. 

5 If one demands that the initial wavelength should equal the island spacing Ip, it 
can be shown that K « Fl% [55]. 
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4-2 Driving force of coarsening 



To gain some insight into the origin of mound coarsening, we consider the 
evolution of the local surface slope m = V/i. The current (56) can be written 
as the negative gradient (with respect to m) of a slope potential defined by 

H 2 

V(m) = -- J duf{u). (63) 
o 



The slope potential allows one to introduce a Lyapunov functional [56] for the 
evolution equation (62), that is, a functional of the surface morphology which 
is a strictly decreasing function of time. Defining 

nh(r,t)} = J d 2 r{^-{V 2 hf + V{Vh)) , (64) 



it is a simple matter to verify that [57,5£ 



* =-y <fr U <o - (65) 



This is an extremely useful relation, because it suggests that the far-from- 
equilibrium growth process can be viewed in analogy to a process of relaxation 
to equilibrium, in which a quantity akin to a free energy is minimized. 

The two terms in the integrand of (64) describe different aspects of this min- 
imization: On the one hand, the value of the slope potential V should be 
minimized locally; on the other hand, the square of the surface curvature V 2 h 
should become small. The minimization of the slope potential drives the pro- 
cess of slope selection. It is evident from (63) that the potential has a minimum 
only if the function / goes through zero at some nonzero slope; otherwise, V 
decreases indefinitely with increasing \m\, and the attempt to minimize it leads 
to unbounded steepening of the morphology. The minimization of the surface 
curvature term in (64) provides, within the continuum equation (62), the driv- 
ing force for coarsening. As will be shown in the next section, this implies that 
the coarsening behavior depends on the spatial distribution of the curvature 
on the surface. 

The analogy of the functional (64) to a thermodynamic free energy becomes 
more pronounced by writing the evolution equation for m in the form 

dm (5^\ . , 

Tsr= vv 'fe)' (66) 
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Fig. 8. Surface morphology generated by numerical solution of the growth equation 
(62) with slope selection (model II). Within the domains the absolute value of 
the slope is |V/i| = 1, and the greyscale encodes the in-plane direction of the height 
gradient. Maxima and minima cannot be told apart, i.e. the morphology is up-down 
symmetric. Courtesy of Martin Rost. 

This is highly reminiscent of the Cahn-Hilliard equation, also known as model 
B, describing the phase ordering of a system with a two-component order 
parameter m [59]. In the case of model II, the order parameter is subject to 
the familiar "Mexican hat" potential which favors \m\ = 1. Mathematically, 
(66) differs from model B in the ordering of the differential operators in front 
of the functional derivative on the right hand side, and in the fact that m 
is irrotational by construction. Physically, an important difference is that the 
domain boundaries between regions with different orientations of the order 
parameter (i.e., the slope) are straight (Figure 8). This is in contrast to 
conventional phase ordering kinetics, where the domain wall curvature drives 
the coarsening process [57,59]. 



4-3 Coarsening laws 

To extract the coarsening behavior from a nonlinear continuum equation such 
as (62), one commonly imposes a scaling hypothesis stating that the patterns 
at different times are similar, in a statistical sense, up to a rescaling by the 
average feature size [59]. In the present context this implies that the height- 
height correlation function (as well as any other statistical measure of the 
morphology) depends on time only through the typical lateral mound size A 



6 More precisely, deformations of the domain boundaries are restored on a time 
scale which is much shorter than the time scale for coarsening [60]. 
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and the surface width W, and hence can be written in the form 

(h(f,t)h(0,t)} = W 2 (t)g(f/X(t)), (67) 

where g is a time-independent scaling function. In addition, the time depen- 
dence of A and W is usually assumed to follow the power laws 

A(t) ~ t n , W ~ t\ (68) 

which define the coarsening exponent n and the roughening exponent (3. The 
mound slope is then of the order of W/ A ~ t@~ n . Slope selection implies (3 = n, 
while in the case of steepening (3 > n. 

To put the scaling hypothesis to work, we consider the time evolution of the 
surface width. Since the mean height has already been subtracted, it is given 
by W 2 = (h 2 ), where the angular brackets represent a spatial average. Multi- 
plying (62) by h and integrating spatially we obtain [61] 

1 dW 2 

^ = (f(\Vh\ 2 )\Vh\ 2 )-((V 2 h) 2 ). (69) 

This clearly demonstrates how the surface is roughened (the width increased) 
by the uphill component of the growth- induced current, and smoothened by 
the curvature term in (62). Both terms on the right hand side of (69) have 
definite signs. Since the coarsening process can be viewed as a competition 
between the two terms, we expect them both, as well as their difference, to 
be of a similar order of magnitude^. This assumption is sufficient to fix the 
exponents n and (3 entering the scaling laws (68). 

According to the scaling hypothesis (67), the typical curvature of the surface is 
of the order of W/X 2 , and hence the right hand side of (69) can be estimated 
as W 2 /X 4 . This immediately leads to n = 1/4 independent of (3. However, 
the estimate implicitly assumes that the curvature is uniformly distributed 
over the surface, which is true for model I (Eq.(57)) without slope selection, 
but not for model II (Eq.(58)). In the presence of slope selection, numerical 
integration of the evolution equation (62) shows that the surface breaks up 
into flat facets at the selected slope, \m\ = 1, which are separated by straight 
domain boundaries [58] (Figure 8). The width of these boundaries introduces 

7 Exact calculations for one-dimensional growth equations show that this assump- 
tion may fail, because the two terms on the right hand side of (69) cancel almost 
completely [62,63]. The simple scaling arguments presented here then only provide 
upper bounds on n and (3. This kind of behavior seems however to be specific to 
the one-dimensional geometry. 
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another length scale into the problem, which invalidates the naive application 
of the scaling hypothesis [59]. This length scale is independent of time, and 
coincides in order of magnitude with the initial mound size 2nV2/q c . The 
surface curvature vanishes on the facets and is concentrated in the network 
of domain boundaries. The correct estimate of the right hand side of (69) is 
therefore ((V 2 /i) 2 ) ~ 1/A, which yields the scaling relation 2(3 — 1 = —n. Since 
slope selection also implies (3 = n, the result n — 1/3 follows for model II. 

To complete the derivation for model I, we note that because of the unbounded 
steepening fi « |V/i|~ 2 for long times, and hence the first term on the right 
hand side of (69) is of order unity. Thus W 2 increases linearly in time, and 
(3 = 1/2. In summary, we have 



The scaling arguments for model I can be extended to equations with a 
smoothening term of the general form — (— l) k (V 2 ) k h with k > 2 [64]. Such 
terms have been proposed to model situations in which thermal detachment 
from steps is impossible, although no clear identification of the associated mi- 
croscopic processes has been provided [41]. The second term on the right hand 
side of (69) then becomes ((V fc /i) 2 ), which is estimated as W 2 /\ 2k , leading to 
n = l/2k. The estimate of the first term remains the same as above, so that 
P — 1/2 independent of k. 

4-4 Crystal anisotropy 

Since the zeros of J determine the stable slopes of the fully-developed mounds, 
the expression for the current should also incorporate the crystal symmetry of 
the surface [50]. To give an example, a possible choice for a surface of square 
symmetry reads [65] 

■h = m x (l -m 2 x - bm 2 ) 



n = 1/4, 13=1/2 (model I) 
n = p = l/3 (model II). 



(70) 
(71) 



J y = m y (l -m 2 



bm 2 x ), 



(72) 



which generates pyramidal mounds with selected slopes 



— * * 

m = 



(±1,±1) 



(73) 



for —1 < b < 1. This reduces to the isotropic model II for b = 1. 
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The microscopic origin of the anisotropy of J is twofold [66]. First, because 
of its dependence on the step structure through, for example, the density of 
kinks, the effective step edge barrier depends on the orientation of the step in 
the plane. Second, the growth- induced current contains contributions from the 
diffusion of atoms along step edges [67,68], which depend even more strongly 
on the step orientation. An approximate evaluation [66] of these contributions 
shows that it is difficult to make contact with phenomenological expressions 
like (72), because the microscopic analysis predicts - in contrast to (72) - that 
the current remains anisotropic even in the limit of small slopes, m — > 0. In 
mathematical terms, this implies a singularity in the function J(m) at rh — 0, 
which may point to a fundamental problem of the continuum theory at small 
slopes. 

An important consequence of crystal anisotropy is the possible breakdown 
of the fundamental scaling hypothesis (67) [65]. For the model (72) this oc- 
curs due to the existence of metastable checkerboard mound patterns, which 
coarsen only through the motion of dislocations ("roof tops"). The distance 
between the dislocations defines a second length scale, which is much larger 
than the mound size. An analysis of the dislocation dynamics suggests that 
n = 1/4, but for quite different reasons than in the case of model I. In the 
presence of hexagonal crystal anisotropy Eq.(71) remains valid, because the 
network of domain boundaries remains sufficiently disordered for the scaling 
arguments to apply [58]. 

4-5 Up-down symmetry and desorption 

As written, Eq.(62) is symmetric under the transformation h — * —h, and 
correspondingly the morphologies it generates are up-down symmetric, that 
is, mounds and valleys have the same shapes. This is in strong disagreement 
with the mound patterns observed in experiments and Monte Carlo simula- 
tions, which display isolated mounds separated by a connected network of 
crevices. These morphologies can be modeled by supplementing (62) by a 
symmetry-breaking term [41]. In a gradient expansion, the lowest order term 
which accomplishes this is of the form V 2 (V/i) 2 [37]. For a one-dimensional 
evolution equation, it has been shown that such a term does not qualitatively 
change the coarsening behavior [62] . Since much of the analysis in this lecture 
relied on the introduction of the slope potential V defined in (63), which be- 
comes impossible in the presence of a symmetry-breaking term, it is unclear 
if the same is true in two dimensions. In fact, it seems that the effect of the 
symmetry-breaking may be quite substantial: If the ridges are removed from 
the network of domain boundaries and only the crevices remain, it is no longer 
possible to localize most of the surface curvature in the domain boundaries, 
and the scaling arguments developed above have to be modified. 
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Fig. 9. Mound configurations generated by numerical solution of the evolution equa- 
tion (74). Left panel shows a greyscale representation of the height. Note the distinct 
up-down asymmetry [69]. 

The up-down symmetry is also broken when desorption from the surface is 
allowed for. The probability for a deposited atom to redesorb from the surface 
before being captured at a step evidently depends on the step density, and 
hence on the surface slope. The desorption rate is therefore an even function 
of V/i. Adding such a function to the right hand side of (62) fundamentally 
changes the character of the evolution equation, because the slope-dependent 
desorption rate cannot be written as the divergence of a current. In [69] the 
effect of desorption on mound coarsening was studied in the framework of the 
evolution equation 

| = -V.[(l-(V^)(V ft )]- TT ^- 5 -(V^, (74) 



where a > is a dimensionless measure of the desorption rate. The numerical 
integration of (74) shows the emergence of conical mounds separated by a 
network of crevices (Figure 9). The form of the desorption term on the right 
hand side of (74) implies that most desorption occurs from maxima and min- 
ima, where Vh ~ 0. Since the minima form a one-dimensional network, while 
the maxima (the tips of the cones) are point-like objects, the growth rate at 
the minima is smaller than that at the maxima by an amount of the order 
of 1/A, where A is the lateral mound size. The surface width then increases 
according to dW/dt ~ 1/A. Together with the slope selection property of (74) 
this implies f3 = n = 1/2. We conclude that desorption leads to a significant 
speedup of coarsening. 
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Fig. 10. Illustration of noise-induced mound coarsening. 
4-6 Noise-induced mound coarsening 

The theory developed so far has been entirely deterministic. Here we show 
that the most important source of fluctuations, the shot noise in the deposition 
beam, induces an alternative coarsening mechanism which generally competes 
with the curvature-driven coarsening described aboveQ [53,70]. 

The basic idea is illustrated in Figure 10. Consider an array of roughly equal 
sized mounds of area A = A 2 and height H = mX/2, where m is the mound 
slope. During a time t, a number F At of atoms is deposited onto a mound, 
with a statistical fluctuation of ±y/FAt. If there is no mass transport between 
neighboring mounds, this translates into a relative height fluctuation 5H = 
V FAt/A. Coarsening then occurs if, by chance, a mound overgrows a less 
fortunate neighbor. The condition for this to happen is that 5H « H, which 
implies the coarsening law 

A « m^e 1 / 4 . (75) 

Under conditions of slope selection, m = const., the coarsening exponent^ is 
n — 1/4, while in general the exponent relation 

p + n = l/2 (76) 

follows. This expresses a competition between coarsening and roughening (or 
steepening): The larger (3, the smaller n, with the limiting case (3 = 1/2, n = 
corresponding to the case of Poisson growth discussed in Section 3. 



8 Quantitative analysis shows, however, that the noise-induced mechanism is neg- 
ligible in most published experiments. 

9 For a (i-dimensional substrate, n = l/(d + 2). 



30 



5 Growth on stepped surfaces 



The orientation of a vicinal surface is close to (in the vicinity of) a high 
symmetry direction of the crystal lattice. Such a surface therefore consists of 
several lattice spacings wide, high index terraces separated by steps, usually 
of monolayer height. During growth on a vicinal surface, the attachment of 
freshly deposited adatoms to the preexisting steps competes with nucleation of 
islands on the terraces. Nucleation is expected to be negligible if the distance 
I between the preexisting steps is small compared to the island spacing l D [71] 
(see also Sect. 4.1). The surface then maintains its vicinal shape, and growth 
occurs through step propagation or step flow. In the following the conditions 
for step flow will be assumed to holdF^j. 

On a perfectly ordered vicinal surface, as it would appear in thermal equi- 
librium at low temperatures, the steps are straight and equally spaced. Cor- 
respondingly, the morphological instabilities of stepped surfaces are of two 
kinds: Either the individual steps develop a meander, beyond their thermal 
or kinetic roughness, or several steps form step bunches, regions of high step 
density separated by large terraces. The main topic of this lecture is a generic 
step meandering instability in homoepitaxial growth, which was first predicted 
theoretically by Bales and Zangwill [73]. It is caused by the asymmetry be- 
tween ascending and descending steps which the step edge barrier introduces. 
Meandering instabilities which have been attributed to the Bales-Zangwill 
mechanism were identified experimentally on surfaces vicinal to Pt(lll) [49] 
and Cu(100) [74-76]. Step bunching during homoepitaxy has been observed 
on several semiconductor surfaces [77,78], but a simple generic mechanism has 
not been suggested. 

5. 1 Stability of a Step Train 

From a theoretical perspective, the step flow growth mode is attractive be- 
cause it can be described in terms of step motion without the need to treat 
island nucleation. The problem simplifies further if the steps are assumed to 
be straight. Then the propagation speed of the j-th step in a step train can 
be written as the sum of the contributions /_ and /+ from the upper and 
lower terraces, each of which is a function of the corresponding terrace width 

10 Due to the stochastic nature of nucleation, the question about the ultimate sta- 
bility of step flow is somewhat subtle [6]. For the one-dimensional Poisson growth 
model of Sect. 3.1, it has been shown that step flow is always metastable, and esti- 
mates for the time scale at which it breaks down have been derived [72]. A similar 
breakdown has been seen in two-dimensional simulations [49], but the underlying 
mechanism is not clear. 
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(Fig. 11). Denoting the position of the j-th step by Xj, the evolution equations 
then read 

dx ' 

-jjT = f+(Xj+l ~ %j) + f-{ x 3 ~ X 'J^)- ( 77 ) 



X J 



X j+1 



Fig. 11. Schematic of a growing step train. 

A train of equally spaced steps moving at speed v = /+(/) + f-(l) evidently 
satisfies these equations. To probe its stability, we consider a small perturba- 
tion of the form 

Xj (t) =jl + v t + e j (t) (78) 

and linearize (77) in the €j. The solutions of the linearized equations are of 
the form ej(t) ~ exp[i<f)j + er (</>)£], where the growth rate o is given in terms 
of the phase shift by the expression 

a(0) = -(1 - cos «/»)(/;(/) - /:(/)) +i sin (f>(f + (l) + f_(l)). (79) 
Stability requires the real part of a to be negative for all 0, which implies 

|(/ + (0-/-(0)>o. (80) 

Roughly speaking, (80) expresses the fact that a step train is stable if the 
steps are fed primarily from the lower terrace, in the sense that /+ > /_ [3]. 
This is easy to understand intuitively: Under this condition a step trailing a 
particularly wide terrace accelerates, and the uniform step spacing is restored. 
When (80) is violated the step train is unstable towards step bunching. The 
largest growth rate is then attained for = tt, hence step pairs form in the 
initial stage of the instability^^] 

While the above analysis applies generally to growing or sublimating vicinal 
surfaces, we now specialize to a surface growing in the absence of evaporation. 
The straightforward evaluation of the functions f± [46] then shows that a 
growing step train is stable whenever D' < D. A step bunching instability 



11 This need no longer be true if long ranged step-step interactions are taken into 
account. 
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during growth would require a "negative" step edge barrier, in the sense of 
D' > D. Similar arguments show that normal step edge barriers generically 
do cause step bunching during sublimation [3,79]. 

The stabilization of the equidistant step train by the step edge barrier may 
be interpreted in terms of an effective, growth-induced step-step repulsion. 
This repulsion is very efficient, in the sense that the resulting terrace width 
fluctuations can be far smaller than in thermal equilibrium [72,80]. 

5.2 The Bales-Zangwill instability 

Bales and Zangwill made the remarkable observation that the very same mech- 
anism that stabilizes a growing vicinal surface against step bunching also 
makes the steps susceptible to a meander instability [73]. Figure 12 illustrates 
the phenomenon on a qualitative level. To account for the mutual repulsion 
between the steps, it is assumed (and will be confirmed by the quantitative 
analysis) that they meander in phase. The terraces can then be subdivided 
into lots, as indicated by the dotted lines. Each lot receives the same number 
of atoms per unit time, which attach primarily to the corresponding segment 
of the ascending step. Because of the meander, the indented segments of the 
step are longer than the protruding ones. Since both capture the same flux, 
the protrusions propagate faster and the deformation is amplified. 



Fig. 12. Schematic of the Bales-Zangwill mechanism for step meandering 

For the quantitative stability analysis, we use the coordinate system shown in 
Fig. 13. The position of the j-th step is described by a function Q(y,t). The 
step spacing of the unperturbed surface is /. Between the steps the adatom 
density n(f, t) satisfies a diffusion equation, which we employ in the stationary 
form 




A 



step 
motion 



DV 2 n + F = 0. 



(81) 
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This is justified provided the time scale for the motion of a step across a 
terrace, 1/F, is large compared to the diffusion time scale l 2 /D, i.e. if the 
Peclet number [81] 

Pe = Fl 2 /D (82) 

is small compared to unity. Since the step spacing in step flow growth has to 
be small compared to Id, which in turn is small compared to (D/F) 1 ^ 2 because 
of the relation (55), the stationarity condition is always satisfiedP* 2 ]. 




x 



Fig. 13. Geometry of the vicinal surface used in the stability analysis. 

The attachment and detachment of adatoms at the steps is described by the 
boundary conditions [79,81] 

± Dn ■ Vn\ x=c .± = k±(n - n^l^-io)- (83) 

Here k + and fc_ denote the attachment rates to an ascending and a descending 
step, respectively, n is the normal vector of the step, and n eq is the equilibrium 
adatom density at the step. The thermodynamic cost of step deformations 
enters the boundary conditions through the expression 

n eq = < (l + ) , (84) 

where 7 = ■y+cPj/d'd 2 is the step stiffness, related to the orientation-dependent 
step free energy 7(1?), n^J is the equilibrium adatom density at a straight step, 
and K st is the step curvature. Equation (84) is the two-dimensional analog of 
the Gibbs-Thomson relation mentioned in Sect. 4.1. 

Once the boundary value problem defined by (81,83,84) has been solved for a 
given configuration of steps, the local normal velocity v%' of each step can be 



Step motion beyond the stationary approximation is treated in [81]. 
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computed from the total mass flux reaching the step from the two adjoining 
terraces, as well as through diffusion along the step. This yields 

= D[n- Vn\ + - n • Vra|_] + J^st^7>4?, (85) 



where /z st is the mobility for migration along a curved (that is, atomically 
rough) step [52], and s denotes the arc length of the step. 

The calculation now proceeds, in principle, as in Sect. 5.1. The general form of 
the perturbed step train is 

( j (y,t)=jl + v t + e j (y,t). (86) 



To linear order in the e^, the solution of the coupled equations can be de- 
composed into normal modes of the form €j(y, t) ~ exp[i(pj + iqy + cr(<f>, q)t\. 
Here q denotes the wavenumber of the step deformation, corresponding to a 
meander wavelength 2ir/q. The real part of the growth rate a(<f),q) turns out 
to be maximal for the in-phase mode = [82]. This is a consequence of the 
kinetically induced step repulsion described in Sect. 5.1: The in-phase mode is 
a compromise which allows the deformed steps to keep the terrace width as 
uniform as possible. 

Since the incipient morphology will be dominated by the fastest growing mode, 
we may assume = in the following. Moreover, for the present purposes it 
is sufficient to consider long wavelength deformations, with a meander wave- 
length large compared to the mean step spacing /. In this limit the expression 
for the growth rate then reads [83] 

Fl*f s 2 (DnWl \ 4 
o-{0,q) = ^—q-[-^^r + ^ s t\iq- (87) 



The positive term proportional to q 2 describes the destabilization of the straight 
step by the attachment asymmetry. The strength of the destabilization is pro- 
portional to the flux F, and to the factor 

fs = {k + k_e/D + k_ + k + ) (88) 



which is a dimensionless measure of the strength of the step edge barrier. The 
negative term proportional to q A describes the thermal relaxation of the step 
towards the (straight) equilibrium shape. The smoothening is driven by the 
step stiffness 7, and it operates through two kinetic channels [84]: Detachment- 
reattachment processes over the terrace, with a rate proportional to the terrace 
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diffusion coefficient and the terrace width, and step edge diffusion with a rate 
proportional to /i st . 



The form of the growth rate (87) is the same that was derived in Sect. 4.1 
for the early stages of the mound instability (Eq.(61)), and correspondingly 
the physics is very similar. For sufficiently small q the quadratic term in (87) 
wins over the quartic term, and hence the step is subject to a long wavelength 
instability for arbitrarily small fluxp 3 "! The range of unstable wavelengths is 
bounded below by 27r/g c , where q c is the wavenumber at which the two terms 
on the right hand side of (87) balance. The dominant meander wavelength Abz 
corresponds to the maximum of (87), which yields Abz = ^V^/Qc- Explicitly, 



\bz 



47T< 



{Dn^l/k B T + /i st )7 
\ FPfs • (89) 



5. 3 Nonlinear step meandering 



An analytic approach to the evolution of the meander instability beyond 
the linear regime has been developed by Misbah, Pierre-Louis and coworkers 
[83,85]. Since the in-phase mode is the most unstable according to linear sta- 
bility analysis, the two-dimensional surface morphology can be represented by 
a one-dimensional function ((y,t) describing the displacement of the common 
step profile from the flat straight reference configuration £ = 0. A solvability 
condition arising from a multiscale expansion in Pe 1 ^ 2 then yields the evolution 
equation 



c* = - 



+ 







0' 



yy 



(90) 



Here subscripts denote derivatives, and the values of the coefficients a, and 
0' can be read off by comparison with (87). The expression in square brack- 
ets is the step curvature, and the terms proportional to and 0' describe 
step smoothening through attachment / detachment kinetics and step edge dif- 
fusion, respectively. The form of these terms follows fro m simp le geometric 
considerations [86]. The two terms differ by a factor of yd + Q, because the 
attachment/detachment kinetics depends on the step width, while step edge 
diffusion does not (compare to (87)). 

Two types of analytic solutions to (90) have been found [85,86]. Station- 
ary solutions are obtained by setting the mass current along the step (the 

13 In the presence of desorption, which is the case originally considered by Bales and 
Zangwill [73], the instability sets in only above a critical flux. 
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quantity inside to curly brackets on the right hand side of (90)) to zero. In 
terms of m(y) = ( y / yjl + Q the stationarity condition reduces to the os- 
cillatory motion of a classical particle in a potential, which can be solved 
by quadratures. One thus obtains a one-parameter family of periodic pro- 
files (s(y) which are most conveniently parametrized by the maximum slope 
S = max y (y. We discuss separately the special cases (3' — and (3 = 0. 
For (3' = 0, the amplitude A(S) of the profile is an increasing function of 
S, while the wavelength A.(S) decreases with increasing S, starting out at 

A(0) = A c = 2n/q c . For S — > oo finite limiting values A(oo) = \j8/3~Ja, 

A(oo) = ^2n(3/a r(3/4)/r(5/4) « 0.5393527.. A c are approached. In contrast, 
for (3 = the potential is harmonic, hence the wavelength of the stationary 
profile is A c independent of its amplitude. 

The separable solutions of interest read [45,85] 

C(i/,t) = 2v / ^err 1 (l-4|t/|/A), -A/2 < y < A/2, (91) 



where the wavelength A is arbitrary. Equation (91) solves (90) exactly in the 
limit t — > oo, when the smoothening terms on the right hand side becomes 
negligible compared to the first term, and the evolution equation reduces to 
(t = —( a /(y)y, the one-dimensional version of the wedding cake equation 
discussed in Sect. 3. 3. The solution (91) is highly singular near the maxima 
and minima, where it diverges as ( ~ ±^/ln(l/|y — ~y~o\), yo = 0, ±A/2. 

40.0 



30.0 
C 20.0 



10.0 



0.0 

0.0 10.0 20.0 30.0 

y 

Fig. 14. Evolution of the step profile starting from a flat initial condition with small 
random fluctuations, for the case of pure step edge diffusion (/3 = 0). Subsequent 
profiles have been shifted in the C-direction. The y-axis has been scaled by A c . 
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In Figure 14 we show results of a numerical integration of (90), starting from a 
small amplitude random initial condition. A regular meander pattern of wave- 
length Abz develops, with an amplitude growing indefinitely as \ft. Closer 
inspection reveals that the sides of the profile follow the separable solution; 
however, the singular spikes at the maxima and minima of (91) are replaced 
by smooth caps which consist of pieces of the stationary solutions. Since the 
slope of (91) increases monotonically upon approaching an extremum while it 
decreases for the stationary profiles, the matching of the two solutions occurs 
near the point of maximum slope. For t — > oo the slope of the separable solu- 
tion diverges, and therefore the cap profile approaches the limiting stationary 
solution Coo(y), an d the length of the cap becomes A(oo). The rescaled step 
profile ((y,t)/\/i approaches an invariant shape in which the cap appears as 
a flat facet. 

5.4 Competing instability mechanisms 

The first quantitative experimental test of the predictions of Bales and Zang- 
will was carried out by Ernst and collaborators for surfaces vicinal to Cu(100) 
[75,76]. Based on the experimentally determined dependence of the meander 
wavelength on temperature and flux, they concluded that their results were 
not consistent with the prediction (89). In particular, the observed flux depen- 
dence A mcan dcr ~ F~ - 2 differs considerably from the predicted Abz ~ F -1 / 2 



This finding suggests that a mechanism different from the one described by 
Bales and Zangwill may be responsible for the step meandering on Cu(100). 
A plausible alternative was proposed in [67], where it was pointed out that a 
one-dimensional analog of the mounding instability discussed in Sections 3 and 
4 should occur on a step, if the diffusion of step adatoms across "descending" 
kinks were suppressed by additional energy barriers. If such barriers are suf- 
ficiently strong, a one-dimensional wedding cake morphology should develop 
along the step, with a characteristic length scale given by the distance 
between the one-dimensional nuclei forming on the straight step in the initial 
stages of growth. This length scale can be estimated from one-dimensional nu- 
cleation theory (see Sect. 2. 3). In order of magnitude, ~ (D e / 'Fid) 1 / 4 , where 
D e is the coefficient of one-dimensional diffusion along a straight (rather than 
kinked) step edge, and F ld = Fl is the effective flux impinging onto a unit 
length of the step edge. A more precise calculation yields [87] 



In contrast to (89), this expression is consistent with the experimental obser- 



[76]. 




(92) 
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vations. 



1000 




Fig. 15. Meander wavelength as a function of flux for a simple cubic solid-on-solid 
growth model. Diamonds and circles refer to conditions with and without facilitated 
edge diffusion, respectively. The full and dashed lines are the predictions of (89) and 
(92) for the two sets of parameters, while the dotted line shows the Bales-Zangwill 
prediction for the case of facilitated edge diffusion. For details on how the material 
parameters entering (89) and (92) are evaluated for the solid-on-solid model, see 
[89]. 

Detailed kinetic Monte Carlo simulations of stepped Cu(100) surfaces support 
the idea of a destabilization of the steps by kink barriers, and show wed- 
ding cake-like step profiles [88]. Simulations of a simple cubic SOS-model have 
moreover demonstrated that both instability mechanisms, the Bales-Zangwill 
scenario and the kink barrier scenario, can be observed in the same system by 
tuning the energy barrier for the detachment of atoms from steps [89]. When 
detachment is suppressed compared to edge diffusion, the kink barrier mecha- 
nism prevails, while the Bales-Zangwill instability is realized when detachment 
rates are comparable to edge diffusion rates. In both cases quantitative agree- 
ment with the predictions (89) and (92) can be achieved without adjustable 
parameters (Figure 15). 
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